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Using perturbation theory in the Euclidean (imaginary time) formalism as well as the nonper- 
turbative Fujikawa method, we verify that the chiral anomaly equation remains unaffected in the 
presence of nonzero chemical potential, fi. We extend our considerations to fermions with exact chi- 
ral symmetry on the lattice and discuss the consequences for the recent Bloch-Wettig proposal for 
the Dirac operator at finite chemical potential. We propose a new simpler method of incorporating 
H and compare it with the Bloch-Wettig idea. 
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Q ■ I. INTRODUCTION 

\^ ' As we know from Noether's theorem, invariance of a Lagrangian of a classical field theory under a continuous 
T-H , symmetry leads to conserved currents. Inclusion of quantum loop corrections can, however, make some currents 
anomalous, and thus lead to the breaking of the corresponding symmetry. Chiral anomalies are a well-known example 
\ of this phenomenon. Chiral anomalies arise in a theory of massless fermions interacting with the gauge fields . 
^ • The flavorless axial current of the fermions is classically conserved but is violated at one-loop level, as was shown 
I ' in the famous calculation of the Adler-Bell-Jackiw(ABJ) triangle diagram for the U{1) case [J, The anomalous 
contribution is a universal feature of the theory and is independent of the ultraviolet regulator used for the quantum 
(-H ■ theory. Fujikawa provided a new insight on anomalies by showing that they arise due to the change of the fermion 
I— I , measure under the corresponding transformation of the fermion fields j5j in the path integral method. Chiral anomalies 
. ■ have a deeper physical significance, as they relate the exact zero modes of the Dirac operator to the nontrivial 
] topological sectors of the gauge fields. Consequently, the chiral anomaly in Quantum Chromodynamics(QCD) is 
QQ ■ thought to give rise to rj' mass 4]. For the physically interesting case of two massless fiavor QCD {Nf — 2), the order 
QQ [ of the chiral phase transition depends Q on the size of the coefficient of the chiral anomaly term. It is of second 
T-H • order, with critical exponents of the 0(4) spin model, if the anomaly contribution is sizeable at finite temperature. 
One could expect a QCD-critical point in the T — fi plane for light quarks in that case. In view of this, it is important 
to ascertain what change occurs in the anomaly in the presence of finite temperature and densities. 
■ In this paper we address both the perturbative and nonperturbative aspects of the chiral anomaly at finite 

^\ , temperature/density. In Sec. I, wc compute the triangle anomaly in the imaginary time formalism of thermal field 
• theory. This method has the advantage that it can be linked to the weak coupling lattice calculations. Lattice QCD 
^ i deals with the imaginary time Euclidean propagators, and hence anomaly calculation in the Euclidean space-time 
would be directly relevant for numerical studies. In Sec. II, we extend Fujikawa's analysis to finite density in the 
continuum. We show that the anomaly equation arising due to the change in the measure of the functional integrals 
under chiral transformation of the fermion fields remains the same at nonzero densities as well. We extend these 
considerations in Sec. Ill to the case of fermions with exact chiral invariance on the lattice. We propose a lattice 
Dirac operator with a term linear in the chemical potential /i, i.e., similar to the continuum and also suggest a way 
to get rid of the spurious divergences in the thermodynamic quantities. Its potential to handle higher order terms in 
the Taylor expansion in chemical potential fi in full QCD is commented upon. 
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FIG. 1: The ABJ triangle diagram(left panel) and its crossed counter part(riglit panel). 



II. ANOMALY AT T = AND /i / IN CONTINUUM 

A. Perturbative calculation 

In this section we calculate the expectation value of the gradient of flavor singlet axial vector current of QCD 
perturbatively in the presence of finite fermion density to check how the anomaly equation is aff'ected in the presence 
of a nonzero chemical potential. The lowest order diagrams are the ABJ triangle diagrams shown in Fig. [TJ It is 
well-known that the higher order diagrams do not contribute to the anomaly equation at zero density, neither do other 
diagrams like the square and pentagon diagrams. We therefore compute only the triangle diagrams at finite density. 
Our starting point is the QCD Lagrangian in the Euclidean space with the finite number density term as defined in 
0. In order to maintain consistency with the lattice literature, we have however chosen the Dirac gamma matrices 
to be Hermitian: 



£ = -ipiP + m)tp - -Tr FapFap + fJ.ilJ^4ip , 



(1) 



where p — ^,j{d„ — ig A'^Ta) with Ta being the generators of the SU(3) gauge group. The ghost terms are not important 
in such a calculation as these do not directly couple to the fermions. The 75 = 71727374 is also Hermitian in our case. 
The inverse free fermion propagator is seen to acquire a ^ dependence and become m + 11^4] . In order to find out 
whether the chiral current j^^ — ipj^j^ip foi' massless quarks is conserved at finite density in one-loop perturbation 
theory, we compute the quantum mechanical expectation value of the derivative of the chiral current i.e. , 



(2) 



where the expectation value of the time ordered product of the three currents at one-loop level is the axialvector- 
vector-vector (AVV) triangle diagram shown in Fig. [TJ Any deviation of this quantity from its classical value would 
give us the anomaly. Using the Euclidean space Feynman rules, the amplitude of the AW triangle diagram can be 
computed. The crossed diagram is the one with the gluon legs exchanged among the two vector (VV) vertices, and 
it corresponds to the process which is quantum mechanically equally favored. 

Denoting by A^'"^(fci, the total amplitude and contracting it with q\, Eq. ([2]) can be written in the momentum 
space for massless quarks as 
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with the tr (Tr ) denoting trace over color (spin) indices. Combining further the first (second) term of the AVV 
diagram and the second (first) term of the corresponding crossed diagram respectively, we rewrite the contracted 
amplitude in terms of functions fi{p, fci) and /2(p, ^2), 



qxA^P'^ = i-i) tT[T'^T']g^ J [hip ~ k,M) ' hip, ^2) + h{p-k^M)~ h{p, ki)] 

where the function fi{p,ki) is defined as, 



(4) 



fi{p,ki) = Tr 



¥1 - «M74 



{pi - +jp _ ku - i/x)2 + {p- kiY 



((P4 - iiJ-Y + P^)((P4 - fci4 - ifJ-Y + (p - kiY) 



since Tr [7V7V7''] = 0. 



(5) 



/2 can be obtained by substituting ^2 for ki and interchanging the indices p and a in Eq. (5). We will use below 
a common notation / for denoting either in order to sketch the proof further. Although the numerator of Eq. ^ 
has terms up to quadratic order in /i, it should be noted that the /i^ terms are ^ /i^Tr ['y5^4^cr^4^pj ^ ^iaAp ^^-^^ 
therefore vanish. In order to further evaluate the right-hand side of Eq. Q , we note that the integrals are linearly 
divergent and hence must be regulated by introducing a cut-off scale, A. This procedure must be carried out in a 
gauge invariant manner such that the vector currents are conserved. In momentum space this amounts to 



A:ipA^'"^(fci,fc2) = A:2.A^''^(fci,fc2 







(6) 



We follow the usual text book 't'I method to impose these conditions above and compute the anomaly. In order to 
highlight the differences due to the fJ. terms, we sketch below the evaluation of just the relevant part of Eq. Q. 
Expanding the first term and combining it with the second, we rewrite the first two integrals as. 



(2^ 



[/(p-fci,fc2)-/(p,fc2)] = /:u- 



-ki^dj + -ki^ki^d^d.f + Oik') 



(7) 



where the derivatives are in the momentum space. The first term of the above integrand can be written as a surface 
integral using Gauss law. 
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(8) 



where we uses the isotropy condition, K^Kd h.^ = gva/'^- It is clear that the second term of the integrand in Eq. ([7]) 
when similarly integrated leads to the gradient of /(p, ^2) at the Fermi surface of radius A, and therefore vanishes as 
Oi-j). Hence this term, and the higher derivative terms, do not contribute in the limit when the cut-off is taken to 
infinity. The other two terms of Eq. as well as the vector current conservation condition Eq. ([SJ, can be similarly 
shown to be fx independent, leading to the canonical result even for fi ^ : 



(9) 



We have thus shown explicitly that the anomaly equation has no corrections due to nonzero fi or, equivalently, at 
nonzero finite density. It is easy to generalize the same computation to nonzero temperatures. At finite temperature, 
the temporal part of the momentum gets quantized as the well-known Matsubara frequencies : p4 = (2"+i)-^ ^ Corre- 
spondingly, ^ gets replaced by ^ where n — ±1, ±2, ±00. The sum over discrete energy eigenvalues, can 
as usual, be split as a zero temperature contribution along with the finite temperature contributions weighted by the 
Fermi-Dirac distribution functions for the particles and antiparticles. Note that the finite temperature contributions 
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will fall off to zero in the ultraviolet limit because these are regulated by the distribution functions. Thus, 
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(10) 



Such perturbative calculations of the AB J anomaly were reported earlier in the real time formalism at finite tempera- 
ture and at both zero and nonzero [1, fermion densities as well as for finite density in Minkowski space-time [ll| . 
We have shown above that these calculations are possible using the imaginary time formalism as well. An imaginary 
time calculation is useful as this can be generalized to weak coupling calculations in lattice gauge theory. 



B. Nonperturbative calculation 



The chiral anomaly in the path integral formalism can also be looked upon as arising due to the change of the 
measure under chiral transformation of the fermion fieldsQ. In this section, Fujikawa's method of anomaly calculation 
in the path integral formalism, at zero temperature and zero fermion density, is extended to the finite fermion density 
case. But before analyzing the finite density problem, the method for /i = is summarized to point out the differences 
that would arise in the finite density case. The partition function for massless fermions interacting with SU {N) gauge 
theory can be written in Euclidean space as 

Z = y" ViPV^[VA^]e- ^ '^^'^ i>J»i>-SYM ^ J 2)^p^[pA^]e-s (11) 

where Sym = 1/2 / d'^x [Tr Fai3{x)Fai3{x) + l/^(/°AJ^)^] is the free Yang-Mills action with appropriate gauge fixing 
f^A"^^ = 0. The action for the ghost term is included within the gauge field measure and hence denoted within square 
brackets. This is justified since we are interested in the change of the fermion fields under chiral transformations and 
the ghost fields do not interact with the fermions. Under the infinitesimal local chiral transformation of the fermion 
fields, given by 

5ijj{x) — ia{x)^^'ijj{x) and 5ij){x) = ia{x)'4){x)^^ , (12) 

the action changes as S ^ S — i J d'^x a{x)d,jjQ. The measure changes as a result of the transformation of the fermion 
fields. The change of measure is, 

V{i;'Vip' = r>V;pVDet| ^y-''^,'' | - DV^PVe"^'-^'^*'' "(x)Tr75 (13) 

where Tr stands for the trace over the spin and the color space. This trace can be computed using the eigenvectors 
of the operator p, since it is an anti-Hermitian operator. It has purely imaginary eigenvalues and the corresponding 
eigenvectors form a complete orthonormal basis. Splitting the trace computation into two parts, the trace over the 
nonzero eigenvalues can be done easily as follows. Since {75, p} = 0, for every eigenvector (pm with nonzero eigenvalue 
A„i ^ 0, there is a corresponding eigenvector 75 0™ with eigenvalue —Am,. Thus for each finite Am , 4>m — 4'm ± 75 0m 
are eigenvectors of 75 with eigenvalues ±1. Since trace is independent of the basis vectors we can also compute the 
trace of 75 in the (jy^ basis. One obtains zero as the result since there are equal number of (jy^ respectively. For the 
zero eigenmodes, p and 75 commute hence each zero mode has a definite chirality, leading to a -1-1 contribution for 
those with 750„ — (j)n and a -1 for the opposite chirality. Hence the complete evaluation of the trace gets a nonzero 
contribution corresponding to the difference between number of the two chiralities: 

Tr75 = ^ (j}\lz<Pn = n+ - n-. (14) 



Chiral Jacobian in the presence of fi 



The presence of finite chemical potential, /z, in the action can be described as an effective change of the Dirac 
operator from p to p — fij^ = p{^i). Under the chiral transformation given in Eq. (jl2l) the action still remains 
invariant as in the zero density case. This is due to the fact that the /i dependent term of the action anticommutes 
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with 75: {75,A*74} = 0. Under the transformations given in Eq. (|12l) the fermion measure changes again by the same 
Jacobian factor Tr75. The corresponding Tr75 is now evaluated in the space of eigenvectors of the new Dirac operator 
This is because the measure is defined by the complete set of states of the Dirac operator which appears in 
the action. Although has both an anti-Hermitian and a Hermitian term, it is still diagonalizable. Consider an 

eigenvector 0,„ of -P(O) with an eigenvalue Am. Let us define two new vectors, Cm and Vm as follows: 

C™(x, r) = e^-0„(x, r) , u/jx, r) = 0t„(x, T)e-^- . (15) 
It is easy to check that Cm is the eigenvector of P{ijl) with the same (purely imaginary) eigenvalue Am, 

PifJ-)Cm = AmCm , (16) 

and vj^ is the eigenvector of |?(/^)^ with the eigenvalue Am — —Am, 

vIPHp) = -X^^vl. (17) 

Note that the sets of eigenvectors {(} and {v} are in one-to-one correspondence with the complete set {(/)}. Using 
the completeness relation for the latter, 

^ /'0m(x,r)0j„(x,r) d*x = I, (18) 

where I denotes the identity matrix, we note that 

^ j Cm(x, T)r;tjx, r) d^x^Y. f r)e^^e-^-cl,l{^, r) d^x - I . (19) 

7n 771 

Moreover, {C} and {v} satisfy the following normality condition, 

vU^, r)Cm(x, r) d^x = [ 0j„e-'^"e''^0,„ d^x - / 4>l{^, r)0^(x, r) d^x = 1 , (20) 



leading to 

z;;f,(x,r)75Cm(x,T) =03,,e^^^75e^^(/),„ = <?!)t,(x, T)75(/)n,(x, r) , (21) 

Using these eigenvector spaces of p{iJ.), the calculation of Tr75 goes through in the same way as for ^(0) above. Since 
the new operator still anticommutes with 75 i.e {75,_p(/i)} = 0, for each eigenvector Cm with eigenvalue Am there is 
an eigenvector 75 Cm with the eigenvalue —Am- Thus the trace of 75 is zero for all nonzero Am- In the basis of the 
zero modes of p{n), given by Cn and wj, the change in the fermion measure is given as, 

Tr75 = E ^i^sCn = ^e-^^Tse'^^n = n+ - n_. (22) 

n n 

Thus the change in the fermion measure due to the chiral transformations is the same as in the zero density case with 
no additional fi dependent terms. Hence the anomaly is unaffected in the presence of /i. Some remarks on the proof 
may be in order. The definition of the vectors Cm and Vm in Eq. (|15|) assumes that neither fi nor r is infinite. The 
same assumption is also utilized in various steps in Eqs. p^ -(P ^ . Clearly at strictly zero temperature, this is not 
tenable. However, an infinitesimally small temperature suffices for the proof to go through. Moreover, since the result 
is finally /.t- independent, we expect the result to be valid at zero temperature, although our proof is valid only in the 
limit of zero temperature. The scaling of the eigenvectors, including the chiral zero modes, by the exp(±/xr) factors 
can be related to a nonunitary transformation of the fermion fields in the QCD action in the presence of /i, given by 

^' (x, r) = e''^ V(x, r) , V^' (x, r) = ^(x, r)c-^^ , (23) 
which makes the action /i- independent: 

5 = [ d^x ^/j'[p - ^74]V'' = / d'^x V^e"^^ - ^74] e'^^V = / d'^x P V - (24) 
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Note that the fields ■0 and -ip at the same space-time point scale difi'erently in the transformation in Eq. ([23)1 which 
is permissible [l^l in the Euclidean field theory since they are mutually independent fields. Let us also emphasize 
that the transformation in Eq. (|23|) is not unitary and thus not physical. Indeed, it merely relates the actions in two 
different physical situations of zero and nonzero /i. One clearly cannot employ it in the evaluation of the partition 
function due to its nonunitary nature. We have shown above that the transformation suggests how to extend the 
cancellation argument for nonzero eigenvalues of the Dirac operator for /i = to the nonzero fj, case as well and is thus 
useful. Furthermore, since the transformation commutes with both flavor singlet and nonsinglet chiral transformations, 
employing it as a prescription to introduce the chemical potential will necessarily lead to a fi dependent action which 
has the same chiral invariance as for ^ = 0. Whether this way to introduce the chemical potential in any theory is 
the only way to do so without affecting its chiral invariance would be interesting to explore; we conjecture that this 
is the case. 



III. ANOMALY ON THE LATTICE AT FINITE DENSITY 



The above discussion of the anomaly in the continuum suggests a way to introduce the chemical potential on the 
lattice. By preserving the transformation ((23| on the lattice, one may expect to maintain the anomaly to remain fi 
independent on the lattice as well. Let us consider the naive massless fermion action on the lattice. 



74 

2 ''x,y+4 



X — i)—6 , " — Uj(x)—S 



■ (25) 



Replacing the ip and ■0 fields in the above action by 0' and respectively, using the lattice analogue of the transfor- 
mation (j23l) . we indeed obtain a fermionic action on the lattice at finite density. 



S 



-^-Ul(x-4)|<5,,^^,-e^ 



Vfu/(x- 1)^,5 Ui(x)^5 



i'y (26) 



with 04 being the lattice spacing in the temporal direction. Unfortunately, the infamous fermion doubling problem is 
related to the fact that the anomaly on the lattice is canceled exactly for such naive fermions. The "no-go" theorem 
of Nielsen and Ninomiya fl^ states that it is impossible to construct lattice Dirac operators which simultaneously 
satisfy Hermiticity, and locality and have chiral symmetry while being free of the "doublers" . The commonly used 
fermions on the lattice, like the Wilson and the Kogut-Susskind fermions do not have ?/a(1) chiral symmetry, and so 
there is no anomaly to speak of. Nevertheless, we note that a similar transformation for such fermions does lead to 
the action popularly used for nonzero chemical potential Iisj . 

Recently, Neuberger |16i| constructed a fermion operator Dpy , commonly known as the overlap operator, which has 
exact chiral symmetry and satisfies the Ginsparg and Wilson|17| relation. 



{l5,Dov} = D0VI5D0V with Dov = 1 + Ib'^i'-fsDw) ■ 
Here e is the sign function and Diy is the canonical Wilson-Dirac operator with a parameter M, 



4 / . 

Dwix, y) = (4 - M) S.,,y -J2(UK^-'^)- 



(27) 



(28) 



The value of the parameter AI is constrained to lie between and 2 for simulating a one fiavor quark on the lattice. 
The overlap fermion action is invariant under the following chiral transformation, as shown by Luscher [isj . 



Stp ^ aj5{l - ^Dov)ip and S^j = aip{l - ^Dov)j5 



(29) 



At zero temperature and density, the change in the measure cornputed on the lattice due to the Luscher transformations 
was shown to be related to the index of the fermion operator [18l - l20l | , 



Tr [275(1 ~ ^Doy)] = -Tr (75^0.) = n+ - = 2 IndexD,, 



(30) 



where n± are right and left handed fermion zero modes respectively. 
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Bloch and Wettig [2l| proposed a method to incorporate the chemical potential in the overlap operator. It consisted 
of i) multiplying U4 by exp(/j,a4)[exp(-/xa4)] in the in Eq- (1^51) and ii) extending the definition of the sign 
function for the resultant complex matrix. The Doy{fj,) also satisfied the Ginsparg- Wilson relation : 

{l5,Dovi^J■)} ^ Dovi^J.h5Dovi^J■)■ with Dovip) = l + -/5e{■^5Dw{^J■)) ■ (31) 

It should be noted that the resultant action does not have the property of eliminating the /i-dependence by any 
transformation like Eg. ((23)) due to the nonlocal nature of Dov 

As we pointed out'22| earlier though, the action S = y '4'x[Dov{l^)]xy'4^y is not invariant under Luscher's chiral 
transformations of Eq. (P^ . Indeed, its variation is easily seen to be 

x,y 

The chiral symmetry violation is of the 0{a) and hence the symmetry is restored in the continuum limit. One may 
alternatively propose modified chiral transformations, 

Sip = a^5{l - -DovitJ-))^' and = - -L'od(m))75 , (32) 

which will ensure SS — 0. In that case, the anomaly equation — Tr (75-Dou(m)) = 2 IndexiDo-u (/i) is valid [2l| on 
the lattice even in the presence of /i, since the fermion measure changes under these transformations by a Jacobian 
factor Tr [275(1 — l/2_Doi)(/^))]- Note, however, that the relevant zero modes are now those of the Dov{^), and thus 
/i dependent, in contrast to our continuum result of the previous section. 

Furthermore, altering the symmetry transformations as above has undesirable physical consequences, as discussed 
in detail in (23j . Let us briefly outline here the main points. Non-Hermiticity of ^^Doy{iJ,) makes the transformations 
nonunitary. The symmetry transformations should not depend on the intensive thermodynamic quantity fi, which is 
a tunable parameter of the physical system. The symmetry group itself changes with /i, leaving no physical order 
parameter which will characterize the chiral phase transition as a function of /i. In contrast, the chiral symmetry 
group remains the same at nonzero temperature (and zero density), allowing us to infer that vanishing of the chiral 
condensate would correspond to restoration of the symmetry for the vacuum. 



A. A simple proposal 

It is well-known that the overlap fermion operator can be obtained [2^, [25j from the five dimensional domain wall 
fermions in the limit of infinite extent of the fifth dimension. The Bloch- Wettig proposal above was also shown to 
arise (26j in this way. It turns out that the chemical potential, ^ enters in their action then as the Lagrange multiplier 
for the number of fermions on each slice of the fifth dimension. This means that all the unphysical "bulk" modes are 
considered with the same weightage in the partition function as the zero modes which eventually correspond to the 
massless quarks in four dimensions. The subsequent cancellation of the bulk contributions using Pauli-Villars fields 
to single out the contribution of a single chiral fermion thus becomes fi dependent on the lattice. Motivated by this 
physical fact, we propose to introduce the chemical potential only to count the fermion confined to the domain wall. 
Integrating out the fermions in the fifth dimension, one is led to the following action, which one would have written 
down naively to add a number density term : 



Dovii^^xy — {-^ov^xy 



2aA M 



(74 + i)c/i(y)^,.,+4 - (1 - 



(33) 



Here Dov is the same Neuberger-Dirac operator of Eq. ([27| . and /t — /ia4 is the chemical potential in lattice units. 
As usual, the volume of the system is ^ = N^a^ and the temperature is T = \l{NTa4) on a x Nt lattice with 
lattice spacings a and 04 in spatial and temporal directions respectively. The term containing the chemical potential 
is not unique. Improved density operators could be used for faster approach to the continuum limit, e.g., addition of 
three-link terms. We could have chosen fl/ s instead of /i/M as the multiplying factor for the conserved number part. 
All such choices of actions are constrained by the fact that these have the correct continuum limit. However the finite 
lattice spacing errors in each of these operators would be different and we comment below on how they may affect 
the numerical simulations. 

Note that our proposal, too, will break exact chiral invariance at the same 0{a) as the Bloch- Wettig proposal. As 
a result, the anomaly equation on the lattice will get -dependent corrections anyway. A significant difference may 
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be the fact that the change in the measure is /i independent for our proposal, as in the case of the continuum. We 
persist with it in the following, nevertheless, as it is simpler and easier to implement. Principally, this is due to the 
fact that one does not have to compute the sign function of a non-Hermitian matrix, with its inherent ambiguities, 
as in the Bloch-Wettig way of incorporating the chemical potential. The non-Hermitian sign function is numerically 
also more expensive to simulate for the full interacting case, whenever that becomes more practical. 

For noninteracting fermions the [/^ = 1 and the above Neuberger-Dirac operator with the chemical potential term 
can be diagonalized in momentum space in terms of the functions, 

hj = —smapj , /14 = sin(a4P4) , 

ft4 

Y^h^ + hj + hl (34) 
such that Dot) (A) can be written as, 

Dov{p,PA, A) = 1 - j7i— - — -77 [74 cos(a4P4) - i sin(a4P4)] . (35) 

s S GaM 

i—l 

To study thermodynamics of fermions one has to necessarily take antiperiodic boundary conditions along the temporal 
direction. Assuming periodic boundary conditions along the spatial directions we obtain 

apj = , rij = 0, .., (TV - 1), 2 = 1, 2, 3 and 

ap4 = w« = ^^""^^^^^ , " = 0, .., {Nt - 1) , (36) 

where aj„ are the Matsubara frequencies. The operator given by Eq. (j35p can be shown to have correct continuum 
limit. The number density can be calculated at zero temperature by the contour integral method as was discussed 
for the Bloch-Wettig version of the overlap fermions at finite in [22| . The major difference one finds is the expected 
/i/a^-divergence (/i^/a^-divergence) in the number (energy) density in the continuum limit of a — > 0. What is perhaps 
not widely appreciated from such calculations is that the leading term, corresponding to the Stefan-Boltzmann limit, 
also changes by a finite computable part. In the next section, we show through numerical evaluations of the sums, 
how one can deal with these problems. 



M 



i=i 



cos apj) (1 — cos (04^4)) ; s 

0,4 
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B. Numerical Results 



We compute two thermodynamic quantities of relevance to the above discussion as well as to the heavy-ion collision 
experiments: the change in the energy density due to the chemical potential, Ae(/i,T) — s{ii,T) — e(0,T) and the 
quark number susceptibility at zero chemical potential, x(0)- These thermodynamic quantities are computed by 
taking appropriate derivatives of the partition function Z = dct£>o„. 

The quantities computed on the lattice are expected to have a dependence on the lattice. In order to 

eliminate these spurious A^ terms, we follow the same prescription which was used for the energy density computation 
at zero temperature (which diverges as A'' ). We compute these thermodynamic quantities at zero temperature and 
subtract them from the corresponding values computed on the lattice at nonzero temperatures. The zero temperature 
values were computed numerically on a lattice with a very large temporal extent Nt and fixed 04 such that T = 
l/^NTdi) — > 0. The Matsubara frequencies then become continuous and hence could be integrated upon numerically. 

Fig. [2] displays the subtracted results for Ae{^,T) for r = fi/T — fiNr — 0.5 and x(0). The former is displayed in 
units of T"* and has the value 0.127 for r = 0.5 in the continuum limit, while the latter is normalized to the ideal gas 
value (T^/3). The M values are as indicated along the symbol used. The subtraction constants had to be cornputed 
separately for energy density and susceptibility. From a comparison of the plots with the corresponding ones [22l | for 
the Bloch-Wettig case, we find that 
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• there are no leftover effects of divergences after the zero temperature subtraction, 

• there are no oscillations for odd-even values of N^, 

• the M-dependence is much less pronounced, and 

• the scaling towards the continuum value is linear with the possibility of an easier extrapolation. 

We also computed the susceptibility using the Wilson fermions and compared the results with those above. We 
found that for Nt = 6 the cut-off effects of the Wilson operator are about 21% larger than the M = 1.60 overlap 
result shown in the right panel of Fig. The difference reduces to about 3% at Nt = 10. Beyond Nt = 10, 

the approach to the continuum limit is almost identical for both the operators. The Wilson fermions have no chiral 
symmetry even for /i = 0, which may make them less favored for the QCD critical point searches which are pivoted 
around the ^ — Q transition. 

We have also checked that there are no other divergent terms of the form 0(a~") with n > 2 in the number density, 
by calculating the fourth-order susceptibility since odd orders of susceptibilities vanish at /i = 0. At zero chemical 
potential, the fourth-order susceptibility is given by, 

X»'(0) = ^C2^i^S^) (38) 



N^Nt V 5/*^ 



A term 0{a ^) in the number density will show up as a divergence in this susceptibility, and will need a subtraction 
too. From Fig. where we display our resuhs for x^^HO) for M = 1.5, we can conclude that there are indeed no 
divergences to be seen in the large iV^ limit. The normalization in this case is also the expected continuum value. It 
is not identical to the Stefan-Boltzmann value of 27r~^. Using the contour integral method it can be easily shown to 
be xi'*^(0) = 2/71-2(1 l/4\ with the additional factor of 0.25 coming from the term usually cancelled in the usual 
prescriptions [13, [H, [13, Hg] . We have found the convergence to the continuum value to be strongly M dependent 
and unfortunately very slow for all values of Af, as seen in the plot B of Fig. ([3]). Introducing the chemical potential 
by choosing ji/ s as the coefficient of the number density term in Eq. (j33l) . instead of the fi/M we used, achieves a 
milder M dependence and a faster convergence towards the continuum. Perhaps improving the number density term 
can achieve a still faster convergence. 



C. A new proposal for QCD critical point via Taylor expansion 

Inspired by the above experience of dealing with the number density in the linear form, as in Eq. p3p . we make 
a proposal valid for all fermions. Because of the infamous sign/phase problem for the fermion determinant with 
nonzero chemical potential, it has been proposed to look for the QCD critical point [1^ by looking for the radius of 
convergence of the Taylor expansion |29l . |30| in fi of the baryonic susceptibility. Computations have been done up 
to the eighth order so far [H, [3l| . Extending these calculations to higher order is both necessary and desirable to 
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FIG. 3: The variation of the ratio of the fourth order susceptibility and the corresponding continuum value as a function of 
l/iVl for C = 4, M = 1.5 for the A) /i/s and B)/i/M ways of incorporating the chemical potential. 



confirm the results already obtained. Our proposal can permit such an endeavor. We denote M (fi) to be any generic 
lattice fermionic operator with the chemical potential /i : 

Sf = 5]4'(x)M(Ai;a;,y)vI/(y) 

= Y,^ix)Dix,y)^iy) + ^^aY,Nix,y) (39) 

x,y x,y 

Here D can be the staggered, overlap, Wilson-Dirac or any other suitable fermion operator, and N(x, y) is the 
corresponding point-split and gauge invariant number density. Eq. ([55]) provides a concrete example of the above for 
the overlap fermions. Note that any improvements in the fermion operator D or the number density TV are generically 
included as long as the classical continuum limit is the same and /i appears linearly. 

It is easy to see that only the first derivative of M with /i is nonzero. All others are zero. Thus denoting by M 
the first derivative of M with respect to and adding more primes in the superscript for successively higher orders, 

M' ^^N{x,y), and M" = M'" = M"" ... = , (40) 

for our proposal to incorporate /i in contrast to the popular exp(±/i) prescription where all derivatives are nonzero: 

M' = M'"... = J2 Nix, y) and M" = M"" = M""" ... ^ . (41) 

x-.y 

As a consequence, the various nonlinear susceptibility expressions, or equivalently the expressions for Taylor series 
coefficients, are a lot simpler and have a lot fewer terms. For example, let us write down a fourth-order coefficient [by 
combining Eqs. (A.4), (A.7), and (A.8) of [2g] : 



y(4) = I_ 
A. y 



Oiiii + 60ii2 + 40i3 + 3O22 + ©4 ) - 3( Oil + O2 ' 



(42) 



Here the notation Oij...i stands for the product, OtOj ■ ■ - Oi. The expressions for On, n=l,4 for our proposal above 
are 

01 = Tr M-^M', (43) 

02 = -Tr M-^M'M-^M', 

03 2 Tr {Nr^M'f, 
Oi = -6 Tr [M-^M'f, 
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in contrast with those for the usual case given in [29| : 

01 = Tr M-H-I', (44) 

02 = -Tr IVr^M'M-^M' + Tt M'^M", 

03 = 2 Tr {IVr^M'f ~3Tt M^^ M' M^^ M" + Tt IVr^M'", 

Oi = -6 Tr {M^^M')^ + 12 Tr {M~^M'fM-^M" - 3 Tr {M^^M")^ 
- 4 Tr M-'^M'M-'^M'" + Tt Nr^M"". 

The eighth-order term needs 0g, which has 18 terms in the usual case whereas it will simply be Os = 
-5040 Tr {M^'^M'f for our proposal. 

The number of matrix inversions required to compute the higher order susceptibilities is also drastically reduced in 
this way of incorporating the chemical potential. This would save a considerable amount of computer time, as matrix 
inversions are the most time intensive operations. Following Fig. 3 of Ref. [HI, one can see that all computations 
referred to on the leftmost branch of the algorithm tree need to be performed when M has a linear /i dependence. 
Thus for the eighth-order susceptibility computation we need to compute only eight matrix inversions as compared 
to the 20 required there, saving 60% of the computer time. For higher order susceptibilities, the number of matrix 
inversions is reduced by at least half, enabling us to compute even higher orders of the Taylor series of thermodynamic 
quantities and thus constrain the radius of convergence and the estimated location of the critical point better. 

Of course, there is a price to pay, and we hope to demonstrate in the future from our ongoing work that it is not very 
big. All the coefficients that one evaluates this way will have the remnants of the terms which are otherwise eliminated 
by the usual prescriptions [13: El; US HI] ■ Based on our computations in the previous section, we suggest that the 
zero temperature contribution to each of them be subtracted by evaluating them on a symmetric N'^ lattice at the 
same (3 — Q/ as the finite temperature computation on the x Nt lattice. Since the second-order susceptibility 
X^^) has a divergence in the continuum limit, its computations may need higher precision to ensure the absence of 
the cut-off effects but the higher order coefficients have no such difficulties. One will also have to rescale the fourth- 
order susceptibility by a factor of 1.25 in order to use it in the ratio or the root method of estimating the radius of 
convergence. We hope that tenth- or even twelfth-order coefficient may thus be computable. 

IV. CONCLUSIONS 

We have shown perturbatively from the computation of the triangle diagram at zero temperature that the anomaly 
equation does not have any finite density correction terms. We have extended our calculations to the nonperturbative 
case where we have used Fujikawa's method to show that the anomaly relation is unaffected in the presence of a 
finite chemical potential. This has an important implication for the lattice field theory in designing the lattice Dirac 
operator for nonzero ^: It should lead to a /^-independent anomaly relation on the lattice. The recent Bloch-Wettig 
proposal for chiral fermion operators at finite density violates the chiral invariance on the lattice itself. While a 
/Lt-dependent modification of the chiral transformation can restore the chiral invariance, it leads to a /i-dependent 
anomaly relation unlike in the continuum theory. Such a modification has other physical consequences discussed in 
Ref. [il. 

We have proposed a physically more justified way of introducing /i in the overlap Dirac operator. In this method the 
chiral symmetry is explicitly broken as well, but the contribution to the anomaly relation from the measure is likely to 
remain independent, with the lattice corrections to the anomaly relation falling off as a power law in the continuum 
limit. It has the expected /i^/a^-type divergences in the continuum limit. We showed how a simple subtraction 
scheme can take care of them in the free case. We proposed to use the simple linear in /i form for the Taylor series 
expansion technique of locating the QCD critical point. It has the advantage that the number of fermion matrix 
inversions goes down drastically when computing the higher order quark number susceptibilities. The higher order 
susceptibility computations are clearly important to accurately locate the critical point in the T-fiB phase space for 
QCD. Our proposal would save much of the computational effort required for obtaining higher order susceptibilities, 
even for the staggered fermions. 
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